%%%% Solution Algorithm for Aiyagari (1994,QJE) Model for Turkey                 
%%%% Written originally by Andreas Mueller in 2015
%%%% Modified by Orhan Torul & Oguz Oztunali
%%%% October 2017, Bogazici University, Istanbul

%%
%%
%% 0. housekeeping
clc;                                % clear screen
clear;                              % clear memory
rng default;                        % fixes the seed for the random number generator (so as to keep simulations the same)
close all;                          % close open windows
%%
%% 1. parameters and functional forms

% parameters
alpha  = 0.56;                      % capital income share from Penn World Table 9.0 (PWT)
b      = 0;                         % borrowing constraint as in Aiyagari (1994, QJE)
beta   = 89/100;                    % subjective discount factor from PWT
delta  = 5.5/100;                   % depreciation rate of physical capital from PWT
gamma  = 3/2;                       % inverse elasticity of intertemporal substitution from macroeconomics literature
varphi = 2/3;                       % Frisch elasticity of labor supply from macroeconomics literature
N      = 5;                         % number of possible productivity realizations

z1     = 96/100;                    % lowest (L) productivity realization
z2     = 98/100;                    % low-medium (LM) productivity realization
z3     = 100/100;                   % medium (M) productivity realization
z4     = 165/100;                   % medium-high (MH) productivity
z5     = 300/100;                   % high (H) productivity

Z  = [z1,z2,z3,z4,z5]';             % vectorizing the productivity grid
Z=Z/1.8;                            % re-scaling parameter for calibration        
%---------------------------------------------------------------------------------------------------------------------
% Markov chain transition probability matrix for productivity

pi      = zeros(N,N);               % initializing the Markov chain transition prob. matrix of NxN

pi(1,1) = 0.70;                     % transition from state 1 (L) to state 1 (L)
pi(1,2) = 0.20;                     % transition from state 1 (L) to state 2 (LM)
pi(1,3) = 1-pi(1,1)-pi(1,2);        % transition from state 1 (L) to state 3 (M)
pi(1,4) = 0;                        % transition from state 1 (L) to state 4 (MH)
pi(1,5) = 0;                        % transition from state 1 (L) to state 5 (H)

pi(2,1) = 0.30;                     % transition from state 2 (LM) to state 1 (L)
pi(2,2) = 0.40;                     % transition from state 2 (LM) to state 2 (LM)
pi(2,3) = 1-pi(2,1)-pi(2,2);        % transition from state 2 (LM) to state 3 (M)
pi(2,4) = 0;                        % transition from state 2 (LM) to state 4 (MH)
pi(2,5) = 0;                        % transition from state 2 (LM) to state 5 (H)

pi(3,1) = 0.10;                     % transition from state 3 (M) to state 1 (L)
pi(3,2) = 0.20;                     % transition from state 3 (M) to state 2 (LM)
pi(3,4) = 0;                        % transition from state 3 (M) to state 3 (M)
pi(3,3) = 1-pi(3,1)-pi(3,2)-pi(3,4);% transition from state 3 (M) to state 4 (MH)
pi(3,5) = 0;                        % transition from state 3 (M) to state 5 (H)

pi(4,1) = 0;                        % transition from state 4 (MH) to state 1 (L)
pi(4,2) = 0;                        % transition from state 4 (MH) to state 2 (LM)
pi(4,3) = 0;                        % transition from state 4 (MH) to state 3 (M)
pi(4,4) = 0.8;                      % transition from state 4 (MH) to state 4 (MH)
pi(4,5) = 1-pi(4,4)-pi(4,3);        % transition from state 4 (MH) to state 5 (H)

pi(5,1) = 0;                        % transition from state 5 (H) to state 1 (L)
pi(5,2) = 0;                        % transition from state 5 (H) to state 2 (LM)
pi(5,3) = 0;                        % transition from state 5 (H) to state 3 (M)
pi(5,4) = 0.05;                     % transition from state 5 (H) to state 4 (MH)
pi(5,5) = 1-pi(5,4)-pi(5,3);        % transition from state 5 (H) to state 5 (H)

rho11=pi(1,1);                      % rho11=prob. of moving to state 1 cond. on state 1
rho12=pi(1,2);                      % rho12=prob. of moving to state 2 cond. on state 1
rho13=pi(1,3);                      % rho13=prob. of moving to state 3 cond. on state 1
rho14=pi(1,4);                      % rho14=prob. of moving to state 4 cond. on state 1
rho15=pi(1,5);                      % rho15=prob. of moving to state 5 cond. on state 1

rho21=pi(2,1);                      % rho21=prob. of moving to state 1 cond. on state 2
rho22=pi(2,2);                      % rho22=prob. of moving to state 2 cond. on state 2
rho23=pi(2,3);                      % rho23=prob. of moving to state 3 cond. on state 2
rho24=pi(2,4);                      % rho24=prob. of moving to state 4 cond. on state 2
rho25=pi(2,5);                      % rho23=prob. of moving to state 5 cond. on state 2

rho31=pi(3,1);                      % rho31=prob. of moving to state 1 cond. on state 3
rho32=pi(3,2);                      % rho32=prob. of moving to state 2 cond. on state 3
rho33=pi(3,3);                      % rho33=prob. of moving to state 3 cond. on state 3
rho34=pi(3,4);                      % rho34=prob. of moving to state 4 cond. on state 3
rho35=pi(3,5);                      % rho35=prob. of moving to state 5 cond. on state 3

rho41=pi(4,1);                      % rho41=prob. of moving to state 1 cond. on state 4
rho42=pi(4,2);                      % rho42=prob. of moving to state 2 cond. on state 4
rho43=pi(4,3);                      % rho43=prob. of moving to state 3 cond. on state 4
rho44=pi(4,4);                      % rho44=prob. of moving to state 4 cond. on state 4
rho45=pi(4,5);                      % rho45=prob. of moving to state 5 cond. on state 4

rho51=pi(5,1);                      % rho51=prob. of moving to state 1 cond. on state 5
rho52=pi(5,2);                      % rho52=prob. of moving to state 2 cond. on state 5
rho53=pi(5,3);                      % rho53=prob. of moving to state 3 cond. on state 5
rho54=pi(5,4);                      % rho54=prob. of moving to state 4 cond. on state 5
rho55=pi(5,5);                      % rho55=prob. of moving to state 5 cond. on state 5

sharesuper=0.25;                    % share of medium-high (MH) plus high (MH) productive individuals


% vectorizing the probability matrix for the equilibrium solver function
rhovector= [rho11 rho12 rho13 rho14 rho15 rho21 rho22 rho23 rho24 rho25 rho31 rho32 rho33 rho34 rho35 rho41 rho42 rho43 rho44 rho45 rho51 rho52 rho53 rho54 rho55];
%---------------------------------------------------------------------------------------------------------------------

% utility and (inverse) marginal utility functions
% contemporaneos utility is defined as u(c)-v(x) with the functinal forms below

u     = @(c) (c.^(1-gamma))/(1-gamma);          % utility from consumption
v     = @(h) (h.^(1+1/varphi))/(1+1/varphi);    % disutility from labor supply
up    = @(c) c.^(-gamma);                       % marginal utility of consumption
invup = @(x) x.^(-1/gamma);                     % inverse of marginal utility of consumption
vp    = @(h) h.^(1/varphi);                     % marginal disutility of labor supply
invvp = @(x) x.^(varphi);                       % inverse marginal disutility of labor supply 


%%
%% 2. discretization

% setting up asset grid
M  = 250;                                       % number of asset grid points
aMax = 100;                                     % maximum asset level
A  = linspace(b,aMax,M)';                       % equally-spaced asset grid from a_1=b *borrowing constraint) to aMax

% vectorizing the grid in two dimensions
Amat = repmat(A,1,N);                           % values of A change vertically
Zmat = repmat(Z',M,1);                          % values of Z change horizontally

% an in-built Matlab built-in alternative can be [Amat,Zmat] = ndgrid(A,Z);

%%
%% 3. endogenous functions

% optimal labor supply
H  = @(c,z,w) invvp(up(c).*z.*w);

% current consumption level, cp0(anext,ynext) is the guess
C0 = @(cp0,r) invup(beta*(1+r)*up(cp0)*pi');
                
% current asset level, c0 = C0(cp0(anext,ynext))
A0 = @(anext,z,c0,r,w) 1/(1+r)*(c0+anext-H(c0,z,w).*z.*w);

%%
%% 4. solving for the stationary equilibrium

% convergence criterion for consumption iteration
crit = 10^(-5);

% parameters of the simulation
% please note that it takes quite some simulation periods to get to the stationary distribution.
% please choose a sufficiently high T >= 10^(4) once the algorithm is running.
I = 10^(4);                                  % number of individuals
T = 1*(10^(4));                              % number of periods

% please choose the interval where to search for the stationary interest rate
% plase note that the stationary distribution is sensitive to the interst rate. 
% please make use of the theoretical result that the stationary rate is slightly below 1/beta-1

r0  = (1/beta-1)-[10^(-12),10^(-2)];
% alternative range as: r0  = (1/beta-1)-[10^(-4),5*10^(-2)];

% setting up an anonymous function searching for equilibrium interest rate
fprintf('Initializing the solution algorithm for the Aiyagari model... \n');
tic;                                        % start timing of the code
myfun   = @(r) to_aiyagari_stationary(r,crit,I,T,Amat,Zmat,alpha,b,delta,rhovector,sharesuper,varphi,A0,C0,H);
options = optimset('display','iter','TolX',1e-5,'MaxIter',2000);
rstar   = fzero(myfun,r0,options);
fprintf('Solved the Aiyagari model in %f sec. \n',toc);

%---------------------------------------------------------------------------------------------------------------------
% retrieving the simulated asset levels
[r,at,ct,ht,inct,zt,labinct,capinct,totalinct,cpnext] = to_aiyagari_stationary(rstar,crit,I,T,Amat,Zmat,alpha,b,delta,rhovector,sharesuper,varphi,A0,C0,H);
htzt=ht.*zt;
lt=htzt;
%---------------------------------------------------------------------------------------------------------------------
% some further preliminaries
r=rstar;
w=(1-alpha)*(alpha/(rstar+delta))^(alpha/(1-alpha));

Zones=ones(M,N);
Zones(:,1)=Zones(:,1)*Z(1);
Zones(:,2)=Zones(:,2)*Z(2);
Zones(:,3)=Zones(:,3)*Z(3);
Zones(:,4)=Zones(:,4)*Z(4);
Zones(:,5)=Zones(:,5)*Z(5);

%---------------------------------------------------------------------------------------------------------------------

% retrieving decision rules

% endogenous labor supply rule
hpnext=H(cpnext,Zones,w);

% endogenous asset choice rule
apnext=zeros(M,N);

for i=1:M
    for j=1:N        
        apnext(i,j)=(1+r)*A(i)+hpnext(i,j)*Z(j)*w-cpnext(i,j);
    end
end

%---------------------------------------------------------------------------------------------------------------------
% plotting decision rules
% graph codes are commented by default so as to keep the running of the code faster

% consumption decision rule
% figure 
% grid
% plot(A,cpnext(:,1),':',A,cpnext(:,2),'--',A,cpnext(:,3),'-.',A,cpnext(:,4),'-.',A,cpnext(:,5),'-.','LineWidth', 3)
% ylabel('Consumption Level','FontSize',16)
% xlabel('Asset Level','FontSize',16)
% legend('Low Productivity (z_1)','Low-Medium Productivity (z_2)','Medium Productivity (z_3)','Medium-High Productivity (z_4)','High Productivity (z_5)','Location','northwest')
% title('Consumption Decision Rule','FontSize',16)
% set(gca, 'FontSize', 16, 'LineWidth', 1, 'Box', 'on', 'FontName', 'Times New Roman');
% xlim([0 45])
% grid
% 
% hours worked decision rule
% figure 
% grid
% plot(A,hpnext(:,1),':',A,hpnext(:,2),'--',A,hpnext(:,3),'-.',A,hpnext(:,4),'-.',A,hpnext(:,5),'LineWidth', 3)
% ylabel('Labor Supply Level','FontSize',16)
% xlabel('Asset Level','FontSize',16)
% legend('Low Productivity (z_1)','Low-Medium Productivity (z_2)','Medium Productivity (z_3)','Medium-High Productivity (z_4)','High Productivity (z_5)','Location','northwest')
% title('Labor Supply Decision Rule','FontSize',16)
% set(gca, 'FontSize', 16, 'LineWidth', 1, 'Box', 'on', 'FontName', 'Times New Roman');
% xlim([0 45])
% %ylim([0 45])
% grid
% 
% 
% % asset choice decision rule
% figure
% plot(A,apnext(:,1),':',A,apnext(:,2),'--',A,apnext(:,3),'-.',A,apnext(:,4),'-.',A,apnext(:,5),A,A,'-','LineWidth', 3)
% ylabel('Next Period Asset Level','FontSize',16)
% xlabel('Asset Level','FontSize',16)
% legend('Low Productivity (z_1)','Medium Productivity (z_2)','High Productivity (z_3)','High Productivity (z_4)','High Productivity (z_5)', '45\circ line','Location','northwest')
% title('Asset Decision Rule','FontSize',16)
% set(gca, 'FontSize', 16, 'LineWidth', 1, 'Box', 'on', 'FontName', 'Times New Roman');
% xlim([0 45])
% ylim([0 45])
% grid



%% 5. generating and reporting distributional results

% histograms

% figure
% % using the last 101 periods
% [nw,xoutw] = hist(at(:,T-100:T),20);     % choosing 20 bins
% bar(xoutw,nw/sum(nw),'FaceColor', [0.75 0.75 0.75]);                   % relative frequency is n/sum(n)
% title('Wealth Distribution','FontSize', 30,'FontWeight', 'bold')
% xlabel('Asset level', 'FontSize', 30,'FontWeight', 'bold')
% ylabel('Relative Frequency','FontWeight', 'bold', 'FontSize', 16)
% set(gca, 'FontSize', 30, 'LineWidth', 1, 'Box', 'on', 'FontName', 'Times New Roman');
% grid
% 
% figure
% [ntinc,xouttinc] = hist(totalinct(:,T-100:T),20);  
% bar(xouttinc,ntinc/sum(ntinc),'FaceColor', [0.45 0.45 0.45]);                  
% title('Income Distribution','FontSize', 16,'FontWeight', 'bold')
% xlabel('Income Level', 'FontSize', 30,'FontWeight', 'bold')
% ylabel('Relative Frequency','FontWeight', 'bold', 'FontSize', 30)
% set(gca, 'FontSize', 30, 'LineWidth', 1, 'Box', 'on', 'FontName', 'Times New Roman');
% grid
% 
% 
% 
% figure
% [nc,xoutc] = hist(ct(:,T-100:T),20);     
% bar(xoutc,nc/sum(nc),'FaceColor', [0.65 0.65 0.65]);                  
% title('Consumption Distribution','FontSize', 16,'FontWeight', 'bold')
% xlabel('Consumption', 'FontSize', 30,'FontWeight', 'bold')
% ylabel('Relative Frequency','FontWeight', 'bold', 'FontSize', 30)
% set(gca, 'FontSize', 30, 'LineWidth', 1, 'Box', 'on', 'FontName', 'Times New Roman');
% grid
% 
% 
% 
% figure
% [nh,xouth] = hist(ht(:,T-100:T),20);   
% bar(xouth,nh/sum(nh),'FaceColor', [0.85 0.85 0.85]);                 
% title('Hours Worked Distribution','FontSize', 16,'FontWeight', 'bold')
% xlabel('Hours Worked', 'FontSize', 30,'FontWeight', 'bold')
% ylabel('Relative Frequency','FontWeight', 'bold', 'FontSize', 30)
% set(gca, 'FontSize', 30, 'LineWidth', 1, 'Box', 'on', 'FontName', 'Times New Roman');
% grid


% please note that in case simulated asset levels take (mildly) negative
% values due to interpolation, one can correct by re-setting them equal to zero. 
% for i=1:10000
%     if asset_distr(i)<=0;
%         asset_distr(i)=0.000001;
%     end
% end


%% calculating stationary distribution values and re-organizing them for reporting
%calculating mean of the simulated series of interest

meanat=mean(mean(at(:,T-100:T)));
meanht=mean(mean(ht(:,T-100:T)));
meanzt=mean(mean(zt(:,T-100:T)));
meanlt=mean(mean(htzt(:,T-100:T)));
meanyt=mean(mean(totalinct(:,T-100:T)));
meanct=mean(mean(ct(:,T-100:T)));
meaninct=mean(mean(inct(:,T-100:T)));
meanlabinct=mean(mean(labinct(:,T-100:T)));
meancapinct=mean(mean(capinct(:,T-100:T)));
meanrt=r;
meanwt=w;

meanmatrix=[meanat meanht meanzt meanlt meanyt+delta*meanat meanct meanrt meanwt meancapinct meanlabinct];

TableStationaryMeans=[meanmatrix];

%exporting stationary equilibrium results to LaTeX
columnlabelStationaryMeans = {'$K$', '$H$' ,'$Z$' ,'$L$', '$Y$', '$C$', '$r$', '$w$' ,'$rK$', '$wL$'};
matrix2latex(TableStationaryMeans, 'TableStationaryMeans.tex', 'columnLabels', columnlabelStationaryMeans,'format', '%-6.3f', 'size', 'footnotesize');


%calculating inequality metrics 
TimeHorizon=100;
GridPointSize=TimeHorizon*I;
distr=ones(GridPointSize,1);        % vectorizing 100 period results
atkinson_epsilon_low=0.25;          % calculating Atkinson index with different epsilon values     
atkinson_epsilon=0.5;
atkinson_epsilon_high=0.75;
atkinson_epsilon_vhigh=1;
atkinson_epsilon_vvhigh=1.5;
atkinson_epsilon_vvvhigh=2;

% wealth distribution metrics
asset_distr=at(:,T-TimeHorizon+1:T);
asset_distr=asset_distr(:);
asset_distr(asset_distr<=0) = 0.0000001;
gini_asset=ginicoeff(asset_distr,distr);
atkinson_asset=atkinsonineq(asset_distr,atkinson_epsilon);
atkinson_asset_low=atkinsonineq(asset_distr,atkinson_epsilon_low);
atkinson_asset_high=atkinsonineq(asset_distr,atkinson_epsilon_high);
atkinson_asset_vhigh=atkinsonineq(asset_distr,atkinson_epsilon_vhigh);
atkinson_asset_vvhigh=atkinsonineq(asset_distr,atkinson_epsilon_vvhigh);
atkinson_asset_vvvhigh=atkinsonineq(asset_distr,atkinson_epsilon_vvvhigh);
theill_asset=theill(asset_distr);
theilt_asset=theilt(asset_distr);
hoover_asset=hoover(asset_distr);

% income distribution metrics
inc_distr=totalinct(:,T-TimeHorizon+1:T);
inc_distr=inc_distr(:);
inc_distr(inc_distr<=0) = 0.0000001;
gini_inc=ginicoeff(inc_distr,distr);
atkinson_inc=atkinsonineq(inc_distr,atkinson_epsilon);
atkinson_inc_low=atkinsonineq(inc_distr,atkinson_epsilon_low);
atkinson_inc_high=atkinsonineq(inc_distr,atkinson_epsilon_high);
atkinson_inc_vhigh=atkinsonineq(inc_distr,atkinson_epsilon_vhigh);
atkinson_inc_vvhigh=atkinsonineq(inc_distr,atkinson_epsilon_vvhigh);
atkinson_inc_vvvhigh=atkinsonineq(inc_distr,atkinson_epsilon_vvvhigh);
theill_inc=theill(inc_distr);
theilt_inc=theilt(inc_distr);
hoover_inc=hoover(inc_distr);

% consumption distribution metrics
c_distr=ct(:,T-TimeHorizon+1:T);
c_distr=c_distr(:);
c_distr(c_distr<=0) = 0.0000001;
gini_c=ginicoeff(c_distr,distr);
atkinson_c=atkinsonineq(c_distr,atkinson_epsilon);
atkinson_c_low=atkinsonineq(c_distr,atkinson_epsilon_low);
atkinson_c_high=atkinsonineq(c_distr,atkinson_epsilon_high);
atkinson_c_vhigh=atkinsonineq(c_distr,atkinson_epsilon_vhigh);
atkinson_c_vvhigh=atkinsonineq(c_distr,atkinson_epsilon_vvhigh);
atkinson_c_vvvhigh=atkinsonineq(c_distr,atkinson_epsilon_vvvhigh);
theill_c=theill(c_distr);
theilt_c=theilt(c_distr);
hoover_c=hoover(c_distr);

% hours worked distribution metrics
h_distr=ht(:,T-TimeHorizon+1:T);
h_distr=h_distr(:);
h_distr(h_distr<=0) = 0.0000001;
gini_h=ginicoeff(h_distr,distr);
atkinson_h=atkinsonineq(h_distr,atkinson_epsilon);
atkinson_h_low=atkinsonineq(h_distr,atkinson_epsilon_low);
atkinson_h_high=atkinsonineq(h_distr,atkinson_epsilon_high);
atkinson_h_vhigh=atkinsonineq(h_distr,atkinson_epsilon_vhigh);
atkinson_h_vvhigh=atkinsonineq(h_distr,atkinson_epsilon_vvhigh);
atkinson_h_vvvhigh=atkinsonineq(h_distr,atkinson_epsilon_vvvhigh);
theill_h=theill(h_distr);
theilt_h=theilt(h_distr);
hoover_h=hoover(h_distr);


ginimatrix=[gini_asset gini_inc gini_c gini_h];
theillmatrix=[theill_asset theill_inc theill_c theill_h];
theiltmatrix=[theilt_asset theilt_inc theilt_c theilt_h];
hoovermatrix=[hoover_asset hoover_inc hoover_c hoover_h];
atkinsonmatrix=[atkinson_asset atkinson_inc atkinson_c atkinson_h];
atkinsonlowmatrix=[atkinson_asset_low atkinson_inc_low atkinson_c_low atkinson_h_low];
atkinsonhighmatrix=[atkinson_asset_high atkinson_inc_high atkinson_c_high atkinson_h_high];
atkinsonvhighmatrix=[atkinson_asset_vhigh atkinson_inc_vhigh atkinson_c_vhigh atkinson_h_vhigh];
atkinsonvvhighmatrix=[atkinson_asset_vvhigh atkinson_inc_vvhigh atkinson_c_vvhigh atkinson_h_vvhigh];
atkinsonvvvhighmatrix=[atkinson_asset_vvvhigh atkinson_inc_vvvhigh atkinson_c_vvvhigh atkinson_h_vvvhigh];

% exporting distribution metrics to LaTeX
TableInequalityMatrix=[ginimatrix; hoovermatrix; theillmatrix; theiltmatrix; atkinsonlowmatrix; atkinsonmatrix; atkinsonhighmatrix; atkinsonvhighmatrix; atkinsonvvhighmatrix; atkinsonvvvhighmatrix];
rowlabelinequality={'Gini Coefficient', 'Hoover Index','Theill Index', 'Theilt Index', 'Atkinson Index $\epsilon=0.25$', 'Atkinson Index $\epsilon=0.50$', 'Atkinson Index $\epsilon=0.75$','Atkinson Index $\epsilon=1.00$','Atkinson Index $\epsilon=1.50$','Atkinson Index $\epsilon=2.00$'};
columnlabelinequality = {'Wealth', 'Income' , 'Consumption', 'Hours Worked'};
matrix2latex(TableInequalityMatrix, 'TableInequalityMatrix.tex','rowLabels', rowlabelinequality, 'columnLabels', columnlabelinequality,'format', '%-6.3f', 'size', 'footnotesize');


%% 6. for welfare cost calculations

% calculating average contemporaneous utilities over individuals & time
for i=1:GridPointSize;
cont_utility_t(i)=u(c_distr(i))-v(h_distr(i));
end
sum_cont_utility=sum(cont_utility_t);
mean_cont_utility=sum_cont_utility/GridPointSize;
